par(mfrow=c(2, 2))
data6_8 <- read.csv("data/习题数据（基于R，EXCEL格式）/csv/E6_8.csv")
ts_x <- ts(data6_8$x, frequency = 12, start=c(1980, 1))

# 时序图
plot(ts_x, type="o", pch=5, col='#39CBB4')



# 1阶12步差分
diff_x <- diff(diff(ts_x, 12))
plot(diff_x)


# ADF检验
# install.packages("aTSA")
library(aTSA)
adf.test(diff_x, nlag=3)

# 白噪声，纯随机检验（<0.05为非白噪声）
for( k in 1:3) print(Box.test(diff_x, lag=6*k, type="Ljung-Box"))

# 自相关图
acf(diff_x, lag.max=30)
pacf(diff_x, lag.max=30)

# 拟合乘法ARIMA模型
fit <- arima(ts_x, order=c(2, 1, 0), seasonal = list(order=c(1, 1, 0), period=12))
fit

# 模型显著性检验
ts.diag(fit)

# 一年预测
# install.packages("forecast")
library(forecast)
x_fore <- forecast(fit, h=12)
x_fore

# 绘制效果图
plot(x_fore, lty=2)
lines(fitted(x_fore), col=4)

